The following files will be used in this lab, all available on Canvas:

  • Climate dataset for Western North America: WNAclimate.csv
  • New York state polygon data in a shapefile: NY_Data.zip
  • A NetCDF file of global monthly air temperature: air.mon.ltm.nc

You will need to download these files from Canvas, and move them from your Downloads folder to the datafiles folder that you made previously. Make sure to unzip the zip files so that R can access the content. Note that on Windows, you will need to right-click on the file and select ‘Extract files’.

You will also need to make sure the following packages are installed on your computer:

pkgs <- c("mapview",
          "raster",
          "RColorBrewer"
          "sf",
          "tmap",
          "viridis")

install.packages(pkgs)
library(mapview)
library(raster)
library(RColorBrewer)
library(sf)
library(tmap)
library(viridis)

Static Maps with tmap

tmap is built on top of ggplot2 and works in a similar way by building a series of layers with map geometries and elements. We start by using tm_shape() to identify the spatial object to be used, and then geometries are added, including filled polygons, borders, legends, etc.


Polygon data

We’ll start by making some maps based on the NY shapefile. Start by loading this, then extract only the polygons belonging to Syracuse.

NY8 <- st_read("../datafiles/NY_data/NY8_utm18.shp")
## Reading layer `NY8_utm18' from data source `C:\Users\u0694980\Box\geog6000\datafiles\NY_data\NY8_utm18.shp' using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## projected CRS:  WGS 84 / UTM zone 18N
Syracuse <- NY8[NY8$AREANAME == "Syracuse city", ]

First, let’s make a simple map showing the polygon outlines using tm_borders():

tm_shape(Syracuse) + tm_borders()

The function tm_fill() can be used to define the polygon fill color based on an attribute in the Syracuse data set (POP8). Note that this automatically adds a legend within the frame of the figure:

tm_shape(Syracuse) + 
  tm_borders() + 
  tm_fill("POP8")

The color scale can be changed by setting the palette argument in tm_fill(). This includes the ColorBrewer scales described above,and the different intervals. For example, to use the ‘Greens’ palette with percentile breaks:

tm_shape(Syracuse) + 
  tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile")

Other map elements can be added. Here we add a longitude/latitude graticule with tm_graticules(), a north arrow with tm_compass(), and a line of text showing the date the map was made with tm_credits().

tm_shape(Syracuse) + 
  tm_graticules(col = "lightgray") + 
  tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile") + 
  tm_compass(position = c("left", "bottom")) + 
  tm_credits("2019-10-19", position = c("right", "top"))

Point data

We’ll next make some maps with the Western North American site data. As before, we load the data, then convert to an sf object. We also load the shapefile of country outlines:

wna_climate <- read.csv("../datafiles/WNAclimate.csv")

wna_climate <- st_as_sf(wna_climate, 
                        coords = c("LONDD", "LATDD"),
                        crs = 4326)

countries <- st_read("../datafiles/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp", 
                     quiet = TRUE)

Individual symbols can be plotted on a color scale using tm_symbols.

tm_shape(wna_climate) + tm_symbols(col = "Jan_Tmp")

This takes the same arguments as tm_fill() for the color palette. We’ll use a red to blue color scale from RColorBrewer. The minus sign before the palette name reverses the order of the colors. As there is a large amount of overlap between the sites, we also add an alpha level to make the symbols transparent.

tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.5, palette = "-RdBu")

We’ll next add country boundaries from the Natural Earth shapefile loaded earlier. Note that as this is a different spatial object, we have to use tm_shape() a second time to reference this, then use tm_borders() to add the lines.

tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") + 
  tm_shape(countries) + 
  tm_borders(col = "gray")

We can also use this

tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") + 
  tm_shape(countries) + 
  tm_borders(col = "gray") + 
  tm_style("cobalt")


Raster data

In the final example, we’ll make figures using the global air temperature dataset. Start by re-reading the data (we also change the name of the layer, to reduce the amount of typing later on…).

air_temp <- rotate(raster("../datafiles/air.mon.ltm.nc", varname = "air"))

air_temp
## class      : RasterLayer 
## dimensions : 73, 144, 10512  (nrow, ncol, ncell)
## resolution : 2.5, 2.5  (x, y)
## extent     : -178.75, 181.25, -91.25, 91.25  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995 
## values     : -37.77505, 33.58953  (min, max)
## z-value    : 0000-12-30
names(air_temp) <- "jan_tmp"

proj4string(air_temp) <- CRS("+init=epsg:4326")

Now let’s make a plot using tm_raster(), which again takes similar options for color palettes. We’ll also add the country borders.

tm_shape(air_temp) + 
  tm_raster(col = "jan_tmp", style = "fisher", palette = "-RdBu") +
  tm_shape(countries) + 
  tm_borders()

We can improve this a little by moving the color legend outside of the plotting area. We’ll increase the number of color classes to 9, and add a histogram showing the frequency of different values.

tm_shape(air_temp) + 
  tm_raster("jan_tmp", 
            style = "fisher", 
            palette = "-RdBu", 
            legend.hist = TRUE, 
            n = 9) +
  tm_shape(countries) + 
  tm_borders() + 
  tm_layout(legend.outside = TRUE, legend.outside.position = "left")


Tips and Tricks

For helpful tips and tricks, try:

tmap_tip()
## An manual legend can be created with tm_add_legend. 
## New since tmap 1.11.2 
## 
## data(metro)
## tm_shape(metro) + 
##   tm_symbols(col = "pop2020", size = "pop2020", legend.col.show = FALSE, legend.size.show = FALSE) + 
##   tm_add_legend(type = "symbol", border.col = "grey40",
##                 col = tmaptools::get_brewer_pal("YlOrRd", 4, plot = FALSE), 
##                 size = sqrt(seq(1:4) / 4),
##                 labels = c("0 to 10 mln", "10 to 20 mln", "20 to 30 mln", "30 to 40 mln"),
##                 title = "Population in 2020")


Interactive maps with mapview

The mapview package provides functions to produce interactive maps in R by building on the R leaflet package, which provides a low-level interface to the javascript library of the same name. mapview supports multiple spatial data types, including points, lines, polygons, and rasters.

mapview(wna_climate, 
        col.regions = "darkred", # fill color
        color = "gray",          # outline color
        alpha.regions = 0.2,     # fill transparency
        alpha = 0.3)             # outline transparency
mapview(Syracuse, 
        col.regions = "darkgreen",
        color = "black")
fn <- system.file("extdata", "kiliNDVI.tif", package = "mapview")

kilimanjaro_NDVI <- raster::stack(fn)[[1]]

mapview(kilimanjaro_NDVI,
        col.regions = viridis::viridis(n = 10)) # supply custom color palette


Here we will only review a few of the many features that mapview provides.

Basemaps

The default basemap called when using mapview() is CartoDB.positron. To use other third-party map tiles, we use the map.types argument in mapview(). For a list of basemap services, see here. For convenience, the leaflet package also provides a named vector of these tile services supported by the leaflet plugin. The vector is called providers and you can access elements of the vector like you would any other:

leaflet::providers[[35]]
## [1] "MapBox"
leaflet::providers$CartoDB.Positron # for the minimalists out there
## [1] "CartoDB.Positron"
mapview(Syracuse, 
        col.regions = "darkgreen",
        color = "black",
        map.types = c("Stamen.Watercolor", "Stamen.TonerLite"))


Layers

mapview provides for layering much in the same way that ggplot2 does.

Syracuse_centers <- st_centroid(Syracuse) 

lyr1 <- mapview(Syracuse, 
                col.regions = "darkgreen",
                color = "black")

lyr2 <- mapview(Syracuse_centers,
                col.regions = "darkblue",
                cex = 3) # size of point

lyr1 + lyr2


Aesthetic Mapping

We can also map attributes of the spatial data to various aesthetics.

mapview(Syracuse,
        zcol = "POP8",
        col.regions = viridis::cividis(n = 10),
        alpha.regions = 0.85)
# breweries in Franconia, Germany
franconia_breweries <- mapview::breweries

# set point size to number of different types of beer served at each brewery
mapview(franconia_breweries, cex = "number.of.types")


Inset Map

For whatever reason, mapview does not currently support interactive inset maps. However, these can still be added using leaflet and a little tinkering with the mapview object.

add_miniMap <- function(x) leaflet::addMiniMap(x@map) 

imap <- mapview(franconia_breweries, cex = "number.of.types")

add_miniMap(imap)